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Abstract 

Maxwell's equations for propagation of electromagnetic waves in dispersive and absorp- 
tive (passive) media are represented in the form of the Schrodinger equation id^>/dt = 
H*$>, where H is a linear differential operator (Hamiltonian) acting on a multi-dimensional 
vector ^ composed of the electromagnetic fields and auxiliary matter fields describing the 
medium response. In this representation, the initial value problem is solved by applying 
the fundamental solution exp(-itH) to the initial field configuration. The Faber polyno- 
mial approximation of the fundamental solution is used to develop a numerical algorithm 
for propagation of broad band wave packets in passive media. The action of the Hamilto- 
nian on the wave function ^ is approximated by the Fourier grid pseudospectral method. 
The algorithm is global in time, meaning that the entire propagation can be carried out 
in just a few time steps. A typical time step is much larger than that in finite differencing 
schemes, Atp ^> . The accuracy and stability of the algorithm is analyzed. The 

Faber propagation method is compared with the Lanczos-Arnoldi propagation method 
with an example of scattering of broad band laser pulses on a periodic grating made of a 
dielectric whose dispersive properties are described by the Rocard-Powels-Debye model. 
The Faber algorithm is shown to be more efficient. The Courant limit for time stepping, 
Ate ~ ll^ll -1 ) is exceeded at least in 3000 times in the Faber propagation scheme. 
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1 Introduction 



Many time-domain algorithms for numerical simulations of the broad band wave packet propa- 
gation in electrodynamics of passive media and/or quantum mechanics use a time stepping, that 
is, given a configuration of the system at time t, a time-domain algorithm produces the system 
configuration at time t + At, where the time step At is determined by conditions resulting from 
the algorithm stability and required accuracy. For instance, in a finite differencing approach, 
such as, e.g., the classical leapfrog scheme, the time step is bounded from above by the stability 
condition (the Courant limit), At < Ate- The upper bound Ate is typically determined by 
the time a signal needs to propagate through an elementary cell of the spatial grid, which is 
by several orders of magnitude smaller than the total propagation time [1]. There is a class 
of problems in numerical electromagnetism where the wave packet dynamics at intermediate 
times is not of significant interest, but rather the final state is important. Computing the scat- 
tering matrix would give one such example. A related and more sophisticated example would 
be simulations of the broad band wave packet propagation in random media [2]. To obtain a 
numerical solution of the initial-value problem in this case, the propagation must be carried 
out multiple times for every (random) state of the medium in order to perform the statistical 
averaging over the medium states. Clearly, a global time-domain algorithm (At ^> Ate) would 
be of great help in reducing computational costs. 

The present work offers a global time-domain algorithm for solving initial value problems 
for Maxwell's equations for passive media whose dispersive and absorptive properties can be 
described by suitable Lorentz, or Rocard-Powels-Debye, or Drude models. The basic idea of 
our approach can be summarized as follows. In Section 2, the Maxwell equations are cast in 
the form of the Schrodinger equation 

("J 

where \& is a multidimensional vector field whose components are electromagnetic fields and 
a set of auxiliary fields that describe the medium response to applied electromagnetic fields 
(e.g., the medium polarization), and H is a linear differential operator that depends on the 
medium dispersive and absorptive properties. Its spectrum is real if no attenuation is present, 
and has a negative imaginary part otherwise. The squared Li norm of \P is proportional to the 
electromagnetic energy of the wave packet. 

If is the initial wave packet configuration, then can be found by using the funda- 
mental solution of Eq. (1.1) 

(t) = e - itH ^ • (1-2) 

Given some (grid) approximation of the spatial dependence of H and \&, Eq. (1.2) provides 
a numerical solution of the initial value problem. In what follows the same letters are used 
for spatial continuum and grid representations of the Hamiltonian and wave functions, unless 
noted otherwise. An exact solution of the initial value problem is understood here in the sense 
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of (1.2) where H is a finite matrix obtained from the continuous Hamiltonian by means of a 
suitable, sufficiently accurate, spatial (grid) representation. 

If H can be directly diagonalized, then (1.2) gives an exact solution for any value of t > 0. 
But this is precisely what one wants to avoid in numerical simulations because the matrix 
H is typically huge and the direct diagonalization is too expensive, if impossible at all. For 
this reason, time domain algorithms use the semigroup property of the fundamental solution: 
exp(-itH) = [exp(— iAtH)] N , where At = t/N with an integer N being the number of time 
steps. For a sufficiently small time step At, typically, At ~ where \\H\\ is the (matrix) 

norm of H, the action of the infinitesimal evolution operator exp(-iAtH) on the state vector 
\1/ can be approximated by various means that do not require any direct diagonalization of H. 

Section 3 is devoted to an algorithm that involves neither a direct diagonalization of H nor 
many time steps. It is based on the well known approximation of an analytical function by 
the Faber polynomial series [3] (see also the textbooks [4]). The Faber approximation method 
has been applied to quantum scattering problems [5] to compute the causal Green's function 
for the Schrodinger equation. The Faber polynomial approximation of the exponential of a 
non-Hermitian operator has also been used to solve the initial value problem for the Liouville 
- von Neumann equation that describes the time evolution of the density matrix in statistical 
systems [6, 7]. In the case when the spectrum of H is real, the approximation yields the well 
known Chebyshev propagation method that has been developed to study wave packet dynamics 
in quantum systems [8, 9, 10] and later used in electrodynamics of non-dispersive media [11]. 

We apply the Faber propagation scheme to solve initial value problems in electrodynamics of 
passive media reformulated in the form of the Schrodinger equation (1.1) with a non-Hermitian 
Hamiltonian, 

n 

tt(f + Atp) = e- iAt * H *(t) » Ck(*t F )F k (H)*(t) . (1.3) 

fc=0 

Here c k (At F ) are the expansion coefficients and F k (H) are Faber polynomials. The action of 
F k (H) on \l/(t) can be computed recursively. The recursion relation depends on the choice of 
the family of Faber polynomials. The latter, in turn, is motivated by spectral properties of H. 
An important point to note is that the expansion (1.3) gives an accurate approximation for the 
fundamental solution for large values of Atp y> Ate ~ 11-^11 1 an d, hence, the propagation can 
be done in just a few time steps. The Faber series (1.3) is known to converge exponentially as 
the approximation order n increases. The accuracy of the algorithm is assessed in Section 4. In 
Sections 5 and 6 the algorithm is applied to scattering of broad band laser pulses on a dielectric 
grating. Dispersive properties of the grating material are described by the Rocard-Powels- 
Debye model with a single pole. The frequency band of the initial pulse is chosen to cover 
the anomalous dispersion range (the pole) of the dielectric. The Faber propagation scheme is 
shown to be more efficient than the Lanczos-Arnoldi propagation scheme applied earlier to the 
same system [12]. The Courant limit can be exceeded in at least 3000 times, Atp > 3000Atc- 
Due to the exponential convergence of the algorithm it can be used as a benchmark for testing 
various time propagation schemes. Note also that it can be applied with any suitable finite- 
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dimensional approximation of the Hamiltonian H (finite elements, or finite differencing, or any 
spectral representation). In our simulations, the Fourier grid pseudospectral representation of 
H has been used [13, 14]. 



2 Maxwell equations in the Hamiltonian form 



Let D and B be electric and magnetic inductions, respectively, and E and H the corresponding 
fields. When no external currents and charges are present, the dynamical Maxwell's equations 
read 

D = cV x H , B = -cV x E . (2.1) 

The over-dot denotes the partial derivative with respect to time, and c is the speed of light in 
the vacuum. Equations (2.1) have to be supplemented by the Gauss law V • D = and also by 
V • B = 0. Relations between the fields and inductions are determined by physical properties 
of the medium in question. 

As an example we consider the Rocard-Powles-Debye model dielectric (the ionic crystal 
model [15, 16]) with one resonance, which is used in our numerical simulations. The case with 
multiple resonances can be studied in a similar fashion. In this model H = B, and the Fourier 
harmonics of the electric field and induction of frequency u are related by D(a>) = e(u)E(u) 
where the dielectric constant is given by 



(so 



£no )UJn 



uj\ — uj 2 — irjuj 



(2.2) 



with Eoofl being constants, ojt the resonant frequency, and i] the attenuation. Let P be the 
dispersive part of the total polarization vector of the medium. Then D = e^E + P. By 
using the Fourier transform, it is straightforward to deduce that P satisfies the second-order 
differential equation 

P + 7]P + c4P = ^oo^pE , 

where up 1 = (e — SooVtAoo ^ £ o ~ * s positive, otherwise, — s 
(2.3) must be solved with zero initial conditions, P = P = 0att = 



(2.3) 

ujp in (2.3). Equation 



Define a set of auxiliary fields Q12 by P = ^/s^uJpQl/^JT and Qi = cj r Q 2 . Maxwell's 
equations and (2.3) can be written as the Schrodinger equation (1.1) in which the wave function 
and the Hamiltonian are defined by 
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(2.4) 



Here e^p are set to one in the vacuum, and to some specific values in the medium in question. 
The squared L 2 norm of the wave function is proportional to the total electromagnetic energy 
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of the wave packet. When attenuation is not present, 77 = 0, the Hamiltonian is Hermitian 
relative to the conventional L 2 scalar product, and the norm (or energy) is conserved. 

In our simulations, an absorbing layer of a conducting medium has been introduced at the 
grid boundaries to prevent reflections of the wave packet. The conductivity a of the layer 
depends on position. The induced current in a conducting media has the form crE. Hence, in 
the presence of the conducting layer the Hamiltonian (2.4) is modified by inserting — Ani^/s^a 
in place of zero in the upper-right corner. Further details can be found in our earlier works 
[17, 18, 12]. 



3 Faber polynomial propagation scheme 

Let D be a bounded, closed continuum in the complex plane such that the complement of D 
is simply connected in the extended complex plane and contains the point at z = 00 (e.g., 
a polygon, an ellipse, etc.). By the Riemann mapping theorem [4], there exists a conformal 
mapping £ which maps the complement of a closed disk with center at the origin and radius p 
onto the complement of D, satisfying the normalization condition, £(w)/w — > 1 as \w\ — > 00. 
Then its Laurent expansion at 00 is given by 

£{w)=w + Y,'YkW- k . (3.1) 

k>0 

The radius p of the disk is called the logarithmic capacity of D. This quantity plays an 
important role in the accuracy analysis given below. The family of Faber polynomials F k 
associated with a conformal mapping £ is defined via the recursion relation 

k 

F k+1 (z) = zF k (z) -Y^ljF k ^(z) - k lk , F (z) = 1 . (3.2) 
3=0 

For a function f(z) that is analytic at every point of D, the Faber series 

00 

f(z)=J2ckF k (z) 

k=0 

is defined by 

1 f /KM) . ,„„> 

where R > p is sufficiently small that / can be extended analytically to the contour Fr being 
the image of the circle \w\ = R under the conformal mapping £. The value R = p is acceptable 
if £ can be extended continuously to the circle \w\ = p (e.g., when the boundary of D is a 
closed simple curve with no self- intersections (a Jourdan curve)). The Faber series converges 
uniformly and absolutely to / on every region bounded by Yr to which / can be extended 
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analytically [19]. This theorem establishes mathematical foundations for the Faber polynomial 
approximation (1.3) of the fundamental solution of (1.1). 

The Faber polynomial algorithm for solving initial value problems for (1.1) is as follows. 
First, choose a (Jourdan) contour T that encloses the spectrum of H. Some criteria for choos- 
ing a contour are discussed in the next section. Second, find the corresponding conformal 
mapping £. In particular, if T is a polygon, this task can be accomplished by the Schwartz- 
Christoffel transformation. For complicated polygons, there is a numerical algorithm to do so 
[20]. Next, the Faber expansion coefficients c^Atp) are computed by means of (3.3) where 
f(z) = exp(— iAtpz). The action of the Faber polynomials of H on *&(£) in (1.3) is computed 
using the recursion relation (3.2). Let = Fk(H)^(t). Then 

k 

= (H - k lk )$ k - lj^k-j , (3-4) 

3=0 

where $o = ^(t) and $i = (H — 7o)$o- The series (1.3) converges uniformly on the entire 
spectral range of H. In order to make the algorithm memory friendly, it is desired to make 
the sequence of 7^ not only finite, but also as short as possible. In Section 5 we apply this 
algorithm to the Hamiltonian (2.4) and choose an ellipse to enclose its spectrum. 

From the numerical point of view, the recursion relation (3.4) is, in general, unstable because 
the minimax norm of Faber polynomials grows rapidly as their order increases, max^i |-Ffc(z)| < 
2p k (see [21]). In other words, the norm of <3>fc would grow exponentially, while the decay of 
Cfc(Ati?) still provides the convergence of (1.3). However, in a numerical implementation of (3.4), 
one might encounter floating point exceptions with a subsequent loss of accuracy. To avoid this 
instability, the Hamiltonian H must be scaled so that its spectrum lies in the domain whose 
logarithmic capacity is one. If f3 is the scaling factor, then exp(— iAtpH) = exp(—iAt s H s ) 
where H s = H//3 and At s = f3Atp. Thus, in the recursion relation (3.4) the scaled Hamiltonian 
H s and the sequence 7^ generated by the conformal mapping (3.1) with p = 1 must be used, 
while the expansion coefficients in the Faber series (1.3) are determined by 

c k (At s ) = -L ^ exp [-iAUie^)} e'^dtp . (3.5) 

Note that the exponential f(z) = exp(—iAt s z) is an analytic function in the entire complex 
plane so that, assuming T to be a Jourdan curve, one can set R = p in (3.3) and use the fact 
that the spectrum of the scaled Hamiltonian lies in a domain with p — 1 and therefore w = e t(fi 
in (3.3). Clearly, the scaling factor f3 must chosen as small as possible to allow for larger time 
steps Atp = At s /(3. 

4 Accuracy and efficiency assessment 

The range Rh of if is a set of complex numbers (^, //"^)/||^|| 2 obtained for all normalizable 
wave functions Here (•, •) denotes a scalar product, and || • || is the norm associated with it. 
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The norm of the resolvent of H is bounded by [22] 

\\{z-H)^\\<[d{z,R H )]-\ (4.1) 

where d(z, z') = \z — z'\ is the distance on the complex plane, and the distance between z and 
a set Rh is defined as mm z / eRH d(z, z'). Let Y be any closed (Jourdan) curve enclosing the 
spectrum of H . Let P n be a polynomial of order n that is used to approximate the fundamental 
solution of (1.1), that is, exp(— itH)^^ ps P n (H)^>Q. By making use of the Cauchy theorem, it 
is straightforward to see that the accuracy of the approximation is bounded by 



-itH 



V -P n (H)y \ 



1 fe- itz -P n {z) , 
1 ^ dz 



J T z — H 

< Cr||*o||max|e- ite -P n (z)| = e n (r)||* || , (4.2) 



where the constant Cr = L T /[2nd{r, Rh)} and L T is the length of T. To find Cr, Eq. (4.1) 
has been used. Note that Cr depends on H and T, but is independent of the approximation 
order n. Hence, it follows from (4.2) that the error of the polynomial approximation of the 
solution of the initial value problem for (1.1) can be made as small as desired because the Faber 
polynomial approximation P n (z) converges to exp(— itz) absolutely and uniformly in D. 

In addition, it is worth noting that the Faber polynomial approximation provides the so 
called "near best" polynomial approximation of an analytic function. By definition, ||/||oo = 
max zeD \f(z)\. The maximum principle for analytic functions states [4] that if / is analytic in 
D and continuous in the closure of D, then |/| cannot attain its maximum at interior points of 
D. According to (4.2) and the maximum principle for functions analytic in D bounded by T, 
the accuracy e n (r) of a polynomial approximation of an analytic function / of a matrix H (in 
our case, f(z) = exp(-itz)) is 

e n (r) = C r ||e- < * 2 -P n (z)|| 00 . (4.3) 

The fundamental theorem for polynomial approximations of functions analytic in the interior of 
D and continuous in D states that there exists a unique best minimax polynomial approximation 
Pi to f, that is, [23] 

ll/-^ / ||oo< HZ-Pnlloo (4-4) 

for any polynomial P n of order n. In practice, it is not easy to find P[. Suppose we choose some 
polynomial approximation, that is, we define a projection operator V n f = P n . In particular, for 
Faber polynomials V n = , and V^f is given by the truncated Faber series. Then it follows 
from the identity f-V n f = f-pf + V n {Pl - f) that 

||/-K/||oo<(l+||^||)||/-P„ / ||oo. (4.5) 

Thus, our polynomial approximation appears to be "near best", provided the norm of the 
projection operator V n is not so large. For Faber polynomials, one can show that [23, 24] 

\\VF\\<-(±\nn + B) (4.6) 



7T \ 7T 2 
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for n > 1. Here B « 1.773 and = J r > 2n and #(2:) is the angle that is made by a 

line tangent to T with the positive real axis. For a convex D, V = 2n by the Radon theorem. 
In our simulations, D is an ellipse, which is convex, therefore 

\\V£\\<9, n< 835. (4.7) 

Equation (4.7) shows that by using the Faber polynomial approximation to / we do not loose 
more than one decimal place in accuracy as compared with the best minimax polynomial 
approximation. In this case, one can also show that [4] 

"-^ £ Ssi«'»i' (4 - 8) 

for any domain bounded by Yr, R > p, to which / can be extended analytically. Thus, the 
Faber series (1.3) converges exponentially as the approximation order n increases. From (4.8) 
some basic principles for choosing the contour Y follow. 

First, because of the exponential convergence of the Faber series, it is desired to make the 
logarithmic capacity p as small as possible. Alternatively, if p is set to one, the scaling factor (3 
must be as small as possible, that is, the contour should enclose the spectrum of H as tight as 
possible. In principle, if the structure of the spectrum of H (or, at least, its range) is roughly 
known, one can find a polygon that tightly encloses the spectrum. The corresponding conformal 
mapping can be computed numerically [20] . The unfortunate feature of this approach is that the 
infinite Laurent series (3.1) is required. Hence, the recursion relation (3.4) becomes memory 
unfriendly in numerical simulations: All the preceding must be kept in the operational 
memory. Thus, when choosing the contour, one should compromise between the approximation 
order and the memory use efficiency of the algorithm [5]. 

Second, if possible, the contour T should not go too far into the upper part of the complex 
plane to avoid the exponential growth of the factor max r | exp(— iAt s z)\ and to allow for larger 
time steps. Note that the necessary accuracy can still be reached, even if the contour goes 
through the upper part of the complex plane, by increasing the approximation order n. The 
latter, however, would lead to a less efficient propagation scheme because more operations per 
time step is required. 



5 The case of an elliptic contour 

Faber polynomials associated with an elliptic contour have the simplest (shortest) recursion 
relation [4]. For this reason this family of the Faber polynomials have been used in many 
aforementioned applications in quantum and statistical mechanics. Here we use the Faber 
polynomials associated with an ellipse to illustrate the Faber propagation scheme in electrody- 
namics of passive media. 
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Consider T being an ellipse (x — Xq) 2 /a 2 + (y — yo) 2 /b 2 = 1 where z = x + iy. The ellipse is 
an image of the circle \w\ — p under the conformal mapping 

Z(w) = w + 7o + 7i/^ , (5-1) 

where a = p + ji/p, b = p — ^/p, and 70 = x Q + iy is the center of the ellipse. The logarithmic 
capacity of an ellipse is p = (a + b)/2 and 7x = p(a — b)/2. We choose p = 1 so that 

71 = 1 - b . 

In this case, the optimization parameters are the scaling factor (5 and the number b + yo that 
determines the factor max r | exp(— iAt s z)\ = exp[At s (b + y )} in the accuracy (4.8) of the Faber 
approximation. 

The recursion relation (3.4) associated with the elliptic contour has only two terms 

$ fc+ i = (H s - 7o )$fc - 7i$ fc _i , k > 1 , (5.2) 
where <3> = and $1 = (H s — 70) $o- The Faber expansion coefficients have the form 

c k (t s ) = (^^j e- iAt ^J k (2t sy /^) , (5.3) 

Here J k is the Bessel function. When computing the integral (3.5) we assume that 71 > (which 
is consistent with the spectral properties of the Hamiltonian H used in our simulations). The 
exponential convergence of the Faber series can easily be seen from the exponential decay of 
the Bessel function for k > 2At sv /77. 

The Hamiltonian (2.4) cannot have eigenvalues with positive imaginary parts, otherwise 
the energy of the wave packet (the squared norm of \l/) would increase with time, which is not 
possible in passive media. Hence, by physical reasons, the spectrum of the Hamiltonian lies 
in the lower half of the complex plane. It is also clear that the spectrum of the Hamiltonian 
is symmetric about the imaginary axis (for every direction in space, there are incoming and 
outgoing waves). Hence, we set xo = 0. The spectrum of H lies in a rectangle [—E m ,E m ] x 
[—v, 0] with E m and v to be determined below. Our strategy is to find an "optimal" ellipse with 
p = 1 that contains a scaled rectangle [—E s , E s ] x [—v s , 0], where E s = E m / f3 and v s = v/(3. 

First, we determine the bounds, E m and v, on the spectral range of H. Let z^ = 
(^,#^)/||^|| 2 be a point in R H . Let H = H — iV where H = (H + W)/2 = H ] and 
V = i(H — H^)/2 = is positive semidefinite. Then 

E m = max Re = max (^, iJ ^)/||^|| 2 = ||#o|| • (5.4) 

Thus, E m is the maximal eigenvalue of H because H is Hermitian. It can be found by the 
standard numerical procedure. If \l/ n = H ^ n -i for n — 1,2, ... for some initial vector \l>o, the 
sequence H^nH/H^-iH converges to the maximal eigenvalue E rn of H as n increases. A rough 
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estimate for E m can also be obtained by noting that the maximal wave vector supported by the 
grid in the Fourier pseudospectral representation is k max = it /a min with a min being the smallest 
grid step (if a non-uniform grid is used). Hence, E m m ck max . Similarly, 

v = max(-Im^) = max V^/H^H 2 = max{ ^n^fe^o max , rj } = ^n^fe^o max , (5.5) 

where a max is the maximal value of conductivity of the absorbing layer. Here we have used 
that fact that V is diagonal and the medium attenuation 77 is small compared to a max . 

By the symmetry, the center of the ellipse is set to coincide with the center of the rectangle, 

7o = -iv 8 /2 . 

An ellipse that contains the rectangle vertices should satisfy the following condition 

a E s 



b Vb 2 - | 7 o| 2 ' 

Since p = 1, a = b — 2. Equation (5.6) relates b and the scaling factor (3. 



(5.6) 



As has been argued above, to increase the time step Atp = At s /(3, the scaling factor f3 
must be minimal. So, one can take b for which (3 attains its minimal value. The smallest (3 is 
reached when 

, / I I \ 2/3 / \ 2/3 

b fho\\ ( v\ ' (5 7) 



a \E S J \2E, 

Observe that if v = 0, that is, if the spectrum of H is real, the optimal ellipse has 6 = and 
(3 = E m . In this case, the Faber polynomial series is nothing but the Chebyshev polynomial 
series. Unfortunately, when v 7^ by making (3 smaller we increase the number b + yo, that is, 
the ellipse gets farther into the upper half of the complex plane and higher orders of the Faber 
approximation are needed to achieve desired accuracy according to (4.8). So, in our simulations 
we take (3 larger than its minimal value and thereby reduce b + yo by making b smaller (see 
next Section for details). 



6 Applications to nanostructured periodic materials 

As an example of possible applications of the present method to photonics, the Faber propa- 
gation scheme associated with an elliptic contour is applied to scattering of broad band wave 
packets on nanostructured periodic materials, the subject of current interest in photonics [25]. 
We consider a grating made of a periodic array of ionic crystal cylinders in vacuum. This 
system has been previously studied by the Lanczos-Arnoldi time propagation scheme [12]. In 
particular, the role of trapped modes (guided wave resonances) and polaritonic excitations in 
transmission and reflection properties of the grating in the infrared range has been elucidated. 
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Apart from illustrating the Faber propagation scheme, our primary interest is to compare its 
efficiency with the efficiency of the Lanczos-Arnoldi propagation scheme. 

The geometry of the system is sketched in the inset of Fig. 2. The system has a translation 
symmetry along one of the Euclidean axes, chosen to be the y axis. It is periodic along the 
x axis with period D g , while the z direction is transverse to the grating. The packing density 
R/Dg = 0.1, where R is the radius of cylinders and D g = 10.8 /an is the grating period. 
The broad band wave packet is represented by a Gaussian pulse that is about 38 fs long and 
has the carrier frequency of 173 meV. It propagates along the z axis and is linearly polarized 
with the electric field oriented along the y axis, i.e., parallel to the cylinders (the so called 
TE polarization). The spectrum of the wave packet is concentrated in a wavelength domain 
A > D g such that the scattering is dominated by the zero diffraction mode (the reflected and 
transmitted beams propagate mainly along the z-axis). A change of variables is used in both x 
(x = fi(xi)) and z (z — ^2(^2)) coordinates to enhance the sampling efficiency in the vicinity 
of medium interfaces so that the boundary conditions are accurately reproduced by the Fourier 
grid pseudospectral method. A typical size of the mesh corresponds to — 17.3D g < z < 15.3D g , 
and — 0.5D g < x < 0.5D g with, respectively, 384 and 64 mesh points. Note that, because of 
the variable change, a uniform mesh in the auxiliary coordinates (xi,x 2 ) corresponds to a non- 
uniform mesh in the physical (x, z) space. The Lanczos-Arnoldi time propagation is carried out 
with a fixed time step Ati = 0.138 fs. The propagation by the Faber method has been done 
with different time steps At F = jAt L , with j = 25, 50, 100, 200, 400, and 1000 (see below). 

The dielectric function of the ionic crystal material is approximated by the single oscillator 
model (2.3). Following the work [15], we chose the parameters representative for the beryllium 
oxide: = 2.99, Eq = 6.6, ujt = 87.0 meV, and the damping r] = 11.51 meV. Thus, for 
D g = 10.8/im two types of resonances can be excited in the system within the frequency domain 
covered by the incident pulse. Structure resonances are characteristic for periodic dielectric 
gratings. They are associated with the existence of guided wave modes [26, 27]. As has been 
demonstrated previously, in the absence of losses, structure resonances lead to 100% reflection 
within a narrow frequency interval(s) for wavelengths A ~ D g . The second type of resonances 
arise because of polaritonic excitations for wavelengths A ~ D T = 2ttc/ut = 26.9 \xm. These 
are associated with substantial energy losses in the ionic crystal material. A detailed discussion 
of the transmission and reflection properties of this grating can be found in [12]. 

Figure la shows the elliptic contour used in our simulations. Its logarithmic capacity is one, 
p — 1, and the corresponding conformal mapping (5.1) reads £(w) = w — 0.005 i + 0.99/iu so 
that b = 0.01. The scaling factor (3 = E m /E s where E m = 0.6468 and E s = 1.7. The shaded 
area is the rectangle [—E S ,E S ] x [— v s ,0] that contains the range of the scaled Hamiltonian 
H s = H//3 as explained in Section 5. The maximum imaginary part of the scaled Hamiltonian, 
v s = 0.01, is consistent with our choice of the absorbing layer. The order n of the Faber 
polynomial approximation is set by the exponential decay of the expansion coefficients (5.3). 
In our simulations, we demand that \ck\ becomes less than 10 -15 for k > n. The behavior of 
|cfc| is shown in Figure lb for time steps Atp = jAtL with j = 50, 200, and 1000. 
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The transmitted signal is collected on the "virtual detector" located at zj, = 3.22 D g behind 
the structure. The zero-order component of the electric field, 



is shown in Fig. 2. 

The existence of a trapped mode (resonance) can easily be inferred from the temporal 
evolution of the electromagnetic field. The main transmitted pulse is clearly visible. It has 
a significant amplitude and a duration about 38 fs. After the main pulse passes the array, 
it leaves behind excited quasistationary modes which loose their energy by radiating almost 
monochromatic waves. By symmetry, the same radiation of quasistationary modes is registered 
in the reflection direction by a detector placed in front of the layer (not shown here). The 
quasistationary mode associated with polaritonic excitations in the ionic crystal has a wave 
length A ~ D F and is short-lived due to the strong absorption of the material at the resonance 
(the anomalous dispersion region). Therefore the observed lasing effect is mainly due to the 
long-lived structure resonance at A ~ D p . The radiation of this mode appears as exponentially 
damped oscillations coming after the main signal. An exponential decay due to a finite lifetime 
of the quasi-stationary mode is clearly seen. The resonance lifetime is in the picosecond range, 
i.e., a thousand times longer than the initial pulse duration. For lossless media, the existence 
of the quasistationary mode(s) leads to a 100% reflection at the resonant frequency, as has 
been discussed in detail in Refs. [12, 28]. Finally, the concept of trapped modes localized 
on successive layers and interacting with each other provides a theoretical framework for light 
propagation in layered structures such as photonic crystal slabs [29]. 

The main results of the paper are summarized in Fig. 3 and Table 1 where we show the 
precision of the Faber propagation scheme and compare its numerical costs with those of the 
Lanczos-Arnoldi scheme. Figure 3 presents a relative error of the time propagation, defined as 
\{E (zd,t) — E re f(zd,t)}/E re f(zd,t)\, where the reference signal E re f(zd,t) is chosen to be the 
result obtained by the Faber propagation scheme with At F — 50 At^. The choice is motivated 
by a higher precision of the Faber scheme (thanks to its exponential convergence) and by the fact 
that the factor maxp | exp(— iAt F z)\ is minimal for the smallest AtF used in our simulations. 
There is no difference between the results with A F = 25At L and A F = 50At L (see below). 

It follows from our results that the Faber propagation scheme has a higher accuracy than 
the Lanczos-Arnoldi propagation at reduced computational costs. The error was saturated at 
10~ 10 value when 10 significant digits in the calculated signal where found to coincide. The 
peaks correspond mainly to the instants of time when the oscillating electric field is close to 
zero. The gain in the propagation efficiency as compared to the Lanczos-Arnoldi scheme is 
twofold. First, a smaller number of actions of H on is needed to obtain ty(t + At F ). In 
the Faber propagation scheme, it is given by the order of the Faber polynomial approximation 
of the fundamental solution, N F = n. In the case of the Lanczos-Arnoldi scheme, the number 
of actions of H on ^>(t) is given by N L = KAt F /AtL where K is the dimension of the Krylov 
space. For the precision shown in Fig. 3, K — 7. Second, as we have already discussed in 
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(6.1) 
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Ref. [12], for a typical size of the mesh as used here, computational costs of acting by H on 
fy(t) are comparable with those of constructing an orthonormal basis for the Krylov space (by 
means of the Arnoldi process) and projecting the Hamiltonian onto the Krylov space (a K x K 
Hessenberg matrix for a non-Hermitian H). This explains an extra factor 2.5 in the fourth 
column of Table 1. For significantly larger sizes of the mesh, in particular, for 3D simulations, 
the computational costs of acting by H on \1> should prevail, and the gain in the computation 
time should simply scale as N L /N F . 

It is worth noting that memory requirements are lower for the present Faber propagation 
owing to the short recursion relation (associated with an elliptic contour). Indeed, in the case 
of the Lanczos-Arnoldi scheme the number of vectors to be kept in the operational memory 
equals K. 

Finally, in the Lanczos-Arnoldi propagation scheme applied to the above system the time 
step AtL exceeds the Courant limit Ate — H-^11 -1 & t least in three times [12]. Therefore, 
the Faber propagation scheme allows one to exceed the Courant limit at least in 3000 times, 
AtF > 3000 Ate as one can see from Table 1. 

Table 1: Numerical costs and efficiency of the Faber polynomial propagation 



Time step Np: number of N^. number of H^J Computation time gain 

units of AtL operations, Faber operations, Lanczos ~ 2.5 N L /N F 



25 


170 


175 


2.5 


50 


290 


350 


3.0 


100 


525 


700 


3.3 


200 


980 


1400 


3.6 


400 


1890 


2800 


3.7 


1000 


4560 


7000 


3.8 



7 Conclusions 

We have shown that the Faber propagation scheme can successfully be used in electrodynamics 
of passive media. The scheme is global in time, that is, it allows one for time steps that 
exceed the Courant limit in a few orders of magnitude. As a point of fact, the propagation can 
actually be carried out in a single time step if the system in question does not have long-lived 
quasistationary modes (as the structure resonance in the example we have considered above). 
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The essential virtue of the scheme is the exponential convergence, which leads to superior 
accuracy as compared to other time domain methods in passive media. The Faber propagation 
scheme can therefore be used as a benchmark, when comparing various propagation methods. 
If the medium is lossless and no absorber is present, the Faber scheme coincides with the 
Chebyshev propagation scheme, whose high accuracy is well known in time domain methods 
in computational quantum physics. 

Another advantage of the present Faber propagation scheme is a relatively low memory 
demand. This, however, is essentially due to an elliptic contour which leads to a family of Faber 
polynomials that are generated by a short recursion relation. For example, the conventional 
leapfrog (time differencing) propagation scheme requires to have two arrays ty(t) and ty(t — At) 
in the operational memory to compute ^(t + At), while in the Faber scheme associated with 
an elliptic contour, a recursive computation of the sum (1.3) requires storing three arrays 
\l/ m (t + Atp), $ m , and $ m _i, where is the series (1.3) with k = 0, l,...,m < n, and 
m — 1, 2, n being the recurrence running index, \l/ m+ i = ^ m + c m+ i$ m+ i. However, the gain 
of the Faber scheme in efficiency and accuracy is enormous. 

It should be noted that we have not explored a further optimization of the present Faber 
propagation scheme because our main goal was to compare it with the Lanczos-Arnoldi prop- 
agation scheme (which was applied to the above system and shown to be more accurate and 
efficient than a typical finite differencing (leapfrog) scheme). For any application, the opti- 
mization should include the following. First, the spread of the spectrum along the real axis 
is essentially determined by the smallest grid spatial step. So, depending on the accuracy de- 
mand, E m can be reduced. Second, the absorbing layer can also be optimized to reduce the 
spread v of the spectrum along the imaginary axis. In addition, one can try to estimate (e.g., 
by perturbation theory) imaginary parts of eigenvalues with large real parts (of order E m ). 
This would lead to a tighter ellipse. Finally, the contour shape itself can also be optimized, 
which, in general, requires a better knowledge of the spectrum of the Hamiltonian. Thus, for 
a specific problem on hands, the Faber propagation scheme can be made even more efficient 
than the simplest example presented in our work. 
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Figure captions 

Fig. la. The elliptic contour used in our simulations. Results are presented on the 
complex plane of the scaled energy (Re E, Im E). The shaded rectangle contains the range of 
the scaled Hamiltonian H s = H//3. Further details are given in the text. 

Fig. lb. The log-log plot of absolute values \ck\ of the expansion coefficients versus k 
for time steps Atp = jAti with j = 50, 200, and 1000 as indicated in the inset of the figure. 

Fig. 2. Electric field of the zero-order transmitted wave as a function of time measured 
in femtoseconds. The signal is registered by a detector placed behind the periodic layer of ionic 
crystal cylinders. The grating geometry is sketched on the inset of the figure. 

Fig. 3. Relative error (defined in the text) for the zero order wave transmitted through 
the periodic layer of ionic crystal cylinders. Results obtained with the Faber propagation scheme 
(full symbols) and the Lanczos-Arnoldi scheme (open circles) are presented as a function of time 
measured in femtoseconds. The time step for the Lanczos-Arnoldi propagation is At L = 0.138 
fs. The time step for the Faber propagation scheme is given by Atp = j AtL, where the 
correspondence between different symbols and the values of j is indicated in the inset of the 
figure. 
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